Libraries:
# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Functions:
# functions -------------------
source("Functions/veganCovEllipse.R")
# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
distance = 'bray',
k = 3,
trymax = 20,
autotransform = FALSE)
## Run 0 stress 0.1272962
## Run 1 stress 0.127296
## ... New best solution
## ... Procrustes: rmse 0.0008255444 max resid 0.002723831
## ... Similar to previous best
## Run 2 stress 0.1272962
## ... Procrustes: rmse 0.0007830489 max resid 0.002595403
## ... Similar to previous best
## Run 3 stress 0.1272961
## ... Procrustes: rmse 0.0007857632 max resid 0.002554142
## ... Similar to previous best
## Run 4 stress 0.1272959
## ... New best solution
## ... Procrustes: rmse 0.0006474293 max resid 0.00211959
## ... Similar to previous best
## Run 5 stress 0.1275742
## ... Procrustes: rmse 0.03456454 max resid 0.1266969
## Run 6 stress 0.1272961
## ... Procrustes: rmse 0.0001239676 max resid 0.0004235327
## ... Similar to previous best
## Run 7 stress 0.1466649
## Run 8 stress 0.1582653
## Run 9 stress 0.1275739
## ... Procrustes: rmse 0.03522942 max resid 0.128812
## Run 10 stress 0.1319007
## Run 11 stress 0.1275737
## ... Procrustes: rmse 0.0347356 max resid 0.1272635
## Run 12 stress 0.1275743
## ... Procrustes: rmse 0.03530544 max resid 0.1290684
## Run 13 stress 0.127296
## ... Procrustes: rmse 0.0001789619 max resid 0.000570497
## ... Similar to previous best
## Run 14 stress 0.127296
## ... Procrustes: rmse 9.379602e-05 max resid 0.0003106422
## ... Similar to previous best
## Run 15 stress 0.1275734
## ... Procrustes: rmse 0.03498567 max resid 0.1280777
## Run 16 stress 0.1272959
## ... New best solution
## ... Procrustes: rmse 7.379182e-05 max resid 0.0002473518
## ... Similar to previous best
## Run 17 stress 0.127296
## ... Procrustes: rmse 0.0001408249 max resid 0.0004555808
## ... Similar to previous best
## Run 18 stress 0.1272962
## ... Procrustes: rmse 0.0002469237 max resid 0.0008214091
## ... Similar to previous best
## Run 19 stress 0.1275735
## ... Procrustes: rmse 0.03493132 max resid 0.1278801
## Run 20 stress 0.1275737
## ... Procrustes: rmse 0.03467607 max resid 0.1270179
## *** Best solution repeated 3 times
# converges, stress of 0.127 - 0.128
set.seed(2)
# 2 dimensions
mali_2d <- metaMDS(rel_mali5,
distance = 'bray',
k = 2,
try = 40,
trymax = 40,
autotransform = FALSE)
## Run 0 stress 0.2013098
## Run 1 stress 0.2421891
## Run 2 stress 0.2164813
## Run 3 stress 0.2248403
## Run 4 stress 0.2013098
## ... New best solution
## ... Procrustes: rmse 7.336217e-05 max resid 0.0002404534
## ... Similar to previous best
## Run 5 stress 0.2166909
## Run 6 stress 0.2013098
## ... Procrustes: rmse 2.20954e-05 max resid 7.187511e-05
## ... Similar to previous best
## Run 7 stress 0.2299679
## Run 8 stress 0.2334747
## Run 9 stress 0.2067192
## Run 10 stress 0.2063475
## Run 11 stress 0.2248403
## Run 12 stress 0.2052525
## Run 13 stress 0.2013098
## ... Procrustes: rmse 8.623799e-05 max resid 0.0002840637
## ... Similar to previous best
## Run 14 stress 0.2248403
## Run 15 stress 0.2013098
## ... Procrustes: rmse 7.687422e-05 max resid 0.0002349284
## ... Similar to previous best
## Run 16 stress 0.2385195
## Run 17 stress 0.2378598
## Run 18 stress 0.2013098
## ... Procrustes: rmse 7.160904e-05 max resid 0.0002354649
## ... Similar to previous best
## Run 19 stress 0.240332
## Run 20 stress 0.271914
## Run 21 stress 0.2063472
## Run 22 stress 0.2013098
## ... Procrustes: rmse 3.826719e-05 max resid 9.76465e-05
## ... Similar to previous best
## Run 23 stress 0.2337596
## Run 24 stress 0.2518641
## Run 25 stress 0.2401184
## Run 26 stress 0.2013098
## ... Procrustes: rmse 6.309312e-06 max resid 1.438008e-05
## ... Similar to previous best
## Run 27 stress 0.2431912
## Run 28 stress 0.2013098
## ... Procrustes: rmse 2.79146e-06 max resid 8.73469e-06
## ... Similar to previous best
## Run 29 stress 0.2076391
## Run 30 stress 0.2013099
## ... Procrustes: rmse 0.0001601721 max resid 0.0005288505
## ... Similar to previous best
## Run 31 stress 0.2013098
## ... Procrustes: rmse 2.72951e-05 max resid 6.009303e-05
## ... Similar to previous best
## Run 32 stress 0.2185542
## Run 33 stress 0.2248403
## Run 34 stress 0.2248403
## Run 35 stress 0.239287
## Run 36 stress 0.2248403
## Run 37 stress 0.2076443
## Run 38 stress 0.2394146
## Run 39 stress 0.2013098
## ... Procrustes: rmse 1.065455e-05 max resid 3.451836e-05
## ... Similar to previous best
## Run 40 stress 0.2013098
## ... Procrustes: rmse 2.107096e-05 max resid 6.965666e-05
## ... Similar to previous best
## *** Best solution repeated 12 times
# stress still too high to represent in 2 dimensions, lowest is ~ 0.2
rm(mali_2d)
ma.veg_enr <- envfit(mali_3d, ma.meta.data_veg)
ma.veg_spc.envr <- envfit(mali_3d, rel_mali5)
# site data -------------------
ma.site.scrs <- as.data.frame(scores(mali_3d, display = "sites"))
ma.site.scrs <- cbind(ma.site.scrs, Treat_Type = ma.meta.data_veg$Treat_Type)
ma.site.scrs <- cbind(ma.site.scrs, Region = ma.meta.data_veg$Region)
# species data -------------------
ma.spp.scrs <- as.data.frame(scores(ma.veg_spc.envr, display = "vectors"))
ma.spp.scrs <- cbind(ma.spp.scrs, Species = rownames(ma.spp.scrs))
ma.spp.scrs <- cbind(ma.spp.scrs, pval = ma.veg_spc.envr$vectors$pvals)
sig.ma_spp.scrs <- subset(ma.spp.scrs, pval<=0.05)
# envr data -------------------
ma.envr.scrs <- as.data.frame(scores(ma.veg_enr, display = "vectors"))
ma.envr.scrs <- cbind(ma.envr.scrs, env.variables = rownames(ma.envr.scrs))
ma.envr.scrs <- cbind(ma.envr.scrs, pval = ma.veg_enr$vectors$pvals)
ma.sig_envr.scrs <- subset(ma.envr.scrs, pval<=0.05)
Really need to figure out how to represent in 2D: figure out the % variation represented by each axis and use the 2 that represent the most variation (can also do multiple 2D graphs)
library(plotly)
fort.mali3 <- fortify(mali_3d) %>%
filter(score == "sites")
plot_ly(x = fort.mali3$NMDS1, y = fort.mali3$NMDS2, z= fort.mali3$NMDS3, type = "scatter3d", mode = "markers", color = ma.meta.data_veg$Treat_Type, symbols = ma.meta.data_veg$Region)
summary(as.factor(ma.meta.data_veg$Treat_Type))
## Control FallRx Harvest MowRx SpringRx
## 5 3 6 3 6
ma.tt <- adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
print(ma.tt)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 4 1.9569 0.31903 2.1082 0.001 ***
## Residual 18 4.1769 0.68097
## Total 22 6.1339 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(pairwiseAdonis)
pairwise.adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
## $parent_call
## [1] "rel_mali5 ~ Treat_Type , strata = Null , permutations 999"
##
## $FallRx_vs_SpringRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.21549 0.10956 0.8613 0.613
## Residual 7 1.75138 0.89044
## Total 8 1.96687 1.00000
##
## $FallRx_vs_MowRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.48859 0.46035 3.4122 0.1
## Residual 4 0.57276 0.53965
## Total 5 1.06134 1.00000
##
## $FallRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.40379 0.20243 1.7767 0.109
## Residual 7 1.59094 0.79757
## Total 8 1.99473 1.00000
##
## $FallRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.39033 0.23712 1.8649 0.118
## Residual 6 1.25581 0.76288
## Total 7 1.64614 1.00000
##
## $SpringRx_vs_MowRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.65697 0.28336 2.7678 0.015 *
## Residual 7 1.66151 0.71664
## Total 8 2.31848 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $SpringRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.4952 0.15598 1.8481 0.061 .
## Residual 10 2.6797 0.84402
## Total 11 3.1749 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $SpringRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.16729 0.0666 0.6422 0.778
## Residual 9 2.34457 0.9334
## Total 10 2.51186 1.0000
##
## $MowRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.72529 0.32578 3.3823 0.014 *
## Residual 7 1.50107 0.67422
## Total 8 2.22636 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $MowRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.84147 0.41918 4.3303 0.014 *
## Residual 6 1.16594 0.58082
## Total 7 2.00741 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Harvest_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.58418 0.21102 2.4072 0.023 *
## Residual 9 2.18412 0.78898
## Total 10 2.76830 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"
#significant: SpringRx vs. MowRx; SpringRx vs. Harvest (borderline at 0.051); MowRx vs. Harvest; MowRx vs Control; Harvest vs Control
library(labdsv)
ma.meta.data_veg$Treat_Group <- NA
ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "FallRx", 1, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "MowRx", 2, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "SpringRx", 3, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Harvest", 4, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Control", 5, ma.meta.data_veg$Treat_Group)
mali_isa <- indval(mali5_avg, ma.meta.data_veg$Treat_Group)
summary(mali_isa)
## cluster indicator_value probability
## QUPR 2 0.8549 0.002
## KAAN 2 0.7534 0.008
## RUSP 2 0.7097 0.005
## SMGL 4 0.6896 0.017
##
## Sum of probabilities = 4.713
##
## Sum of Indicator Values = 6.68
##
## Sum of Significant Indicator Values = 3.01
##
## Number of Significant Indicators = 4
##
## Significant Indicator Distribution
##
## 2 4
## 3 1
gr <- mali_isa$maxcls[mali_isa$pval <= 0.05]
iv <- mali_isa$indcls[mali_isa$pval <= 0.05]
pv <- mali_isa$pval[mali_isa$pval <= 0.05]
fr <- apply(mali5_avg > 0, 2, sum)[mali_isa$pval <= 0.05]
fridg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fridg <- fridg[order(fridg$group, -fridg$indval),]
print(fridg)
## group indval pvalue freq
## QUPR 2 0.8549332 0.002 6
## KAAN 2 0.7534400 0.008 10
## RUSP 2 0.7096551 0.005 7
## SMGL 4 0.6896051 0.017 6
# 1 is FallRx
# 2 is MowRx
# 3 is SpringRx
# 4 is Harvest
# 5 is Control